2019年开发竞赛[ENVI遥感应用组]特等奖作品欣赏——青藏高原冻土监测制图
欢迎大家关注公众号“ENVI技术殿堂”
作品单位：中国地质大学（武汉）地理与信息工程学院
小组成员：陈一仰 范增辉 黄博文 高瑞翔
指导教师：刘修国 张正加
获奖情况：特等奖
https://v.qq.com/txp/iframe/player.html?width=500&height=375&auto=0&vid=y3045n6fgao
一、作品概述
冻土是指零摄氏度以下，并含有冰的各种岩石和土壤。一般可分为短时冻土（数小时/数日以至半月）、季节冻土（半月至数月）以及多年冻土（又称永久冻土，指的是持续二年或二年以上的冻结不融的土层）。
目前，全球气候变化导致冻土退化情况日益严重，冻土的冻结和融化的反复交替会造成地质环境的破坏，从而导致房屋和道路等地面工程建筑物的地基破裂或者塌陷，还会引起山体滑坡、洪水暴发以及冰川移动等问题。
青藏高原是地球上面积最大的高海拔冻土区，地质地貌独特，对气候变化敏感，易受人类活动影响。根据研究表明，20世纪后期青藏高原总的冻土面积退化了约23.84％。由于青藏高原海拔高，地域广环境恶劣等特点，给实地调查带来困难。而遥感手段具有监测范围广，分辨率高，重返周期短等优点，可得到最新的冻土情况并通过变化趋势研究冻土对不同气候因子的响应情况和冻土对环境的影响，对生态研究有着重要意义。
本应用结合多源遥感数据基于TTOP(Top temper of permafrost)模型得到中国青藏高原地区不同尺度的冻土分布图，根据历史冻土分布数据和钻孔数据对制图成果进行评价，并结合地面站点的降雨、气温及海拔数据分析气候因素对冻土变化的影响，最后使用C#+ArcEnging+IDL进行平台开发。
二、作品主要技术流程
本应用整体流程框图划分为四个层次：数据层、处理层、应用层与分析层。数据层集成多源遥感数据，处理层对数据进行预处理、地表覆盖类型分类与温度反演；应用层承接上层结果，进行冻融指数与土壤导热系数计算，基于TTOP模型反演活动层顶板温度进行冻土制图，最后分析层结合历史与实测数据进行精度评价，分析冻土对气候变化的时空响应关系，探究冻土退化趋势。
图1 整体应用流程图
2.1 遥感数据预处理
本作品通过编程开发实现对于遥感数据的批量重投影、裁剪、拼接等预处理。
2.2 基于PSC的温度反演
Landsat数据采用Practical single-channel(PSC)算法进行温度反演，PSC算法是在单窗算法的基础上，针对不同传感器计算固定参数，降低了单窗算法中由于普朗克函数线性化与大气校正的带来的误差。该算法分两步进行LST估计，首先是计算黑体辐射亮温，通过同时消除发射率和大气效应计算地表辐射度(B(Ts))。再通过反演普朗克函数，由B(Ts)计算LST。此种方法为了避免由普朗克函数线性化引起的回归不确定性。
2.3 土壤导热系数计算
通过植被覆盖系数对应表生成土壤导热系数，之后通过对30*30像素范围内的所有土壤导热系数进行加权平均重采样，获得尺度变换后的导热系数。
2.4 冻融指数计算
2.4.1 年极值月温法冻融指数计算
对于Landsat数据冻融指数计算，根据实验区附近站点地温数据，求解年最高温度，最低温度所在月。假定全年温度变化符合以下曲线，每年使用极值温度所在月的影像，根据公式即可得到最终结果。
图2 温度拟合曲线
2.4.2 全年数据累加法冻融指数计算
传统方法假设温度服从简谐波动，根据温度的年最大和最小月平均值计算融化指数与冻结指数。然而，由于全球气候变化的影响，气候波动和极端气候出现的次数有所增加，使得年内温度变化不再符合简单的简谐波动假设，造成积温计算的不准确。针对这一问题，本作品采用一种基于全年温度产品的冻融指数计算方法，以0℃为界限，将一年的温度划分为正温和负温，并分别累加计算正积温和负积温，从而得到准确的冻结指数和融化指数。
2.5 基于TTOP模型的冻土分布反演方法
TTOP模型是用于反演活动层顶层的温度的模型。根据冻融指数和导热系数即可进行相应的计算。计算公式为:
根据此公式，输入土壤导热系数以及冻融指数即可得到活动层顶层的温度。将每9年温度平均处理，根据程国栋提出的热力学分类系统(表1)对青藏高原冻土进行分类。
表1：热力学分类系统
冻土类型 | 年均气温(℃) | TTOP(℃) |
极稳定型 | <-8.5 | <-5 |
稳定型 | -8.5～-6.5 | -5～-3 |
亚稳定型 | -6.5～-5 | -3～-1.5 |
过渡型 | -5～-4 | -1.5～-0.5 |
不稳定型 | -4～-2 | -0.5～0.5 |
极不稳定型 | -2～-1 | >0.5 |
2.6 精度评价
对2001-2009年冻土分布图和2010-2018年冻土分布图分别与1:300万青藏高原冻土分布图数据进行宏观上的叠加对比，比较这两者冻土分布的大致范围，从而在宏观上确认冻土分布的精确性。同时将冻土分布图与地面钻孔数据对比，得出MAGT值与实测钻孔数据差值均值为-0.618，标准差为0.930。表明制图结果能够较好地反映出冻土区域的实际情况。
图3 2001-2009年冻土分布图叠加
图4 2010-2018年冻土分布图叠加
2.7 联合多气候要素的冻土时空变化分析
冻土时间变化分析，对青藏高原2001-2018年CRU（英国东英格利亚大学气候研究所）全球月平均气温与降水数据进行统计，结果发现2001-2018年整体上气温呈缓慢上升的趋势，其中2010-2018年的平均气温要比2001-2009年的平均气温高0.0143℃；降水呈缓慢上升的趋势。其中2010-2018年的平均降水量要比2001-2009年的平均降水量多0.8577mm。
图5 2001-2018青藏平均气温变化图
图6 2001-2018年降水变化图
冻结指数、融化指数代表该地区的气温的冻结能力与融化能力，而融化指数的加强与冻结指数的下降势必会导致冻土的退化。如下图所示，2001-2018年青藏高地区融化指数在缓慢上升，而冻结指数有下降的趋势，且融化指数上升速度要稍大于冻结指数下降的速度。
图7 2001-2018年DDT与DDF变化图
对于冻土空间变化分析，将青藏高原海拔高程数据与2001-2018年每年冻土情况叠加分析，统计每类冻土对应的平均海拔，得出不同类型的冻土在海拔分布上也呈现出了逐渐升高的趋势，各类冻土在平均海拔分布上存在着很强的线性相关关系：
表2 冻土类型与海拔关系
使用冻土类型转移矩阵分析冻土的退化情况，统计两期（2001到2009年为一期，2010到2018为二期）之间的土地利用转移矩阵。可以发现绝大多数的永久冻土在退化且主要集中在一个级别之内，冻土越不稳定变化的可能性越高。根据表4的面积统计也可以发现，除不稳定型冻土之外，其他的冻土面积均在减少，且冻土类型越稳定面积减少的越多。
表3 冻土类型转移矩阵
图8 冻土变化示意图
表4 近18年来青藏高原多年冻土类型面积统计（×104km2）
三、作品成果展示——冻土监测制图平台
结合应用所需要，本作品冻土监测制图开发主要涉及到五大功能模块。分别为通用工具模块、温度反演模块、地表覆盖分类模块、冻土分布模块、制图输出模块。
通用工具模块：实现了多种栅格的数据的读取、矢量数据的读取、以及.xmd工程文件的读取。在视图方面实现了鹰眼视图。同时支持放大、缩小、全图显示以及漫游的功能。并且可以更具用户的喜好更改平台样式。
图9 基础功能模块
温度反演模块：本作品设定了一键温度反演功能，只需要输入Landsat8影像，MODIS水汽数据，试验区矢量数据就可直接生成温度产品。
图10 温度反演模块
冻土反演模块：实现对MODIS温度产品的预处理、计算冻融冻结指数、根据地表分类求取地面导热系数和冻土成图这几部分功能。
图11 冻土分布成图模块
制图输出模块：实现了基本的制图要素的显示，图例、指北针、比列尺、标题。
图12 制图输出模块
根据上述功能操作，本作品可以一键制图获得MODIS温度数据以及青藏高原冻土分布数据。
图13 温度制图
图14 冻土分布制图
四、作品关键技术与亮点
(1)本作品运用Practical single-channel(PSC)算法进行温度反演，在单窗算法的基础上针对不同传感器计算固定参数，降低了单窗算法中由于普朗克函数线性化与大气校正的带来的误差。同时根据全年数据累加法计算冻融指数，避免了传统方法中对温度服从简谐波动的假设，一定程度上消除了气候波动对年内积温计算准确性的影响。
(2)根据历史冻土分布数据与钻孔数据进行精度评价，表明本作品冻土制图结果具有较好精度。
(3)结合CRU全球气温降水数据、高程数据对冻土时空变化进行趋势分析，建立冻土变化转移矩阵，对冻土“由何变，变为何”有了清楚直观的掌握，为青藏高原地区冻土保护提供一定决策支持。
(4)本作品开发的青藏高原冻土监测制图一体化平台，针对青藏高原冻土反演制图中的关键步骤和主要问题，集成了温度反演、冻土反演和制图输出三大模块。温度反演模块实现了基于PSC算法的一键温度反演；冻土反演模块能够根据MODIS地温产品和地表覆盖分类实现冻土反演；制图输出模块能够一键输出温度专题图、MAGT专题图、多年期冻土专题图等。该系统功能完善，算法准确，高度自动化地实现了冻土参数反演制图。